Non-Gaussian Resistance Fluctuations in Disordered Materials 
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We study the distribution of resistance fluctuations of conducting thin fihns with different levels 
of internal disorder. The film is modeled as a resistor network in a steady state determined by the 
competition between two biased processes, breaking and recovery of the elementary resistors. The 
fluctuations of the film resistance are calculated by Monte Carlo simulations which are performed 
under different bias conditions, from the linear regime up to the threshold for electrical breakdown. 
Depending on the value of the external current, on the level of disorder and on the size of the system, 
the distribution of the resistance fiuctuations can exhibit significant deviations from Gaussianity. 
As a general trend, a size dependent, non universal distribution is found for systems with low 
and intermediate disorder. However, for strongly disordered systems, close to the critical point of 
the conductor-insulator transition, the non- Gaussianity persists when the size is increased and the 
distribution of resistance fluctuations is well described by the universal Bramwell-Holdsworth-Pinton 
distribution. 
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I. INTRODUCTION 

In a macroscopic system, the fluctuations of a global 
quantity are generally expected to follow a Gaussian dis- 
tribution around the average value. This expectation 
relies on the assumption, frequently true, that the sys- 
tem can be divided into statistically independent mi- 
croscopic or mesoscopic elements. For such a system, 
the conditions of validity of the central-limit theorem 
are in fact satisfied and, consequently, Gaussian distri- 
butions of fluctuations are actually observed. On the 
other hand, when the correlations between the differ- 
ent elements of the system are significant, the central- 
limit theorem is no longer valid and a non-Gaussian dis- 
tribution of fluctuations is expected. Therefore, non- 
Gaussian fluctuation distributions of global quantities are 
the signature of a correlated system. As correlations be- 
come important near the critical points of phase tran- 
sitions, non-Gaussian distributions are usually observed 
near criticality^^^. In these conditions, the self-similarity 
of the system over all the scales, from a characteristic mi- 
croscopic length up to the size of the system (which for 
a system of finite size sets the cut-off for the correla- 
tion length), has important implications on the fluctua- 
tion distribution^"'*'^'^"'''^. Recently, new light has been 
shed on this subject by the remarkable discovery made by 
Bramwell, Holdsworth and Pinton (BHP)^ of a common 
behavior of the distribution of fluctuations of two quite 
different systems. Namely, the distribution of the fluc- 



tuations of the power consumption measured in confined 
turbulent-flow experiments and the distribution of mag- 
netization fluctuations in a two-dimensional XY model in 
the spin- wave regime at low temperature'^. Since these 
two systems appear to have very little in common, BHP 
advanced the reasonable suggestion that the origin of this 
common distribution should be attributed to scale invari- 
ance, the only property apparently shared by the two 
systems^. This suggestion has been supported by the 
subsequent finding that many scale invariant systems, 
in both equilibrium and non-equilibrium steady states, 
display the same functional form for the distribution of 
fluctuations^'^'^°'^^'^. Concerning this point, it must be 
noted that the probability density function (PDF) of the 
distribution of fluctuations of a critical system can be 
considered as universal, in the sense that, when properly 
normalized, it depends on very few basic symmetries that 
define the universality class of the system*. Thus, the 
BHP distribution describes a universal behavior of the 
fluctuations of a critical system. In any case, it must be 
remarked that not all critical systems fluctuate accord- 
ing to this distribution'^^^'^-'^. On the other hand, far from 
criticality, the correlations between different elements of 
the systems can also be important. This is particularly 
true for systems in non-equilibrium steady-states, where 
non-Gaussian fluctuations are frequently observed^'^'^'^^. 
Therefore the study of non-Gaussian fluctuations and the 
understanding of their link with other features of the sys- 
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tcm can provide new insights on some basic properties of 
complex systems^'^'*"'^'"'^"'^'^. 

Here, we study the distribution of resistance fluctua- 
tions of conducting thin films of different size and with 
different levels of internal disorder. The resistance fluc- 
tuations are investigated under different bias conditions, 
from the linear response regime up to the threshold for 
electrical breakdown. This last phenomenon, which con- 
sists of an irreversible increase of the resistance, is thus 
associated with a conductor-insulator transition and it 
occurs in conducting materials stressed by high current 
dcnsitics"'^^"-'^^. In our study we make use of the Station- 
ary and Biased Resistor Network (SBRN) model^°^^^. 
This model provides a good description of many features 
associated with the electrical instability of composites 
materials^'''^* and with the electromigration damage of 
metal lines^^^^^, two important classes of breakdown phe- 
nomena. The film is modeled as a resistor network which 
reaches a steady state determined by the competition 
between two biased stochastic processes, breaking and 
recovery of the elementary resistors. The resistance and 
its fluctuations are then calculated by Monte Carlo sim- 
ulations. Resistance fluctuations are found to deviate 
from Gaussianity near electrical breakdown. As a gen- 
eral trend, a non universal distribution of fluctuations 
is found for systems with low and intermediate disorder, 
where the deviations from Gaussianity vanish in the large 
size limit. However, for highly disordered systems, close 
to the critical point of the conductor-insulator transition, 
these deviations from Gaussianity persist when the size is 
increased and the distribution of resistance fluctuations 
is well described by the universal BHP distribution^''*. 



II. MODEL 

We describe a conducting film with granular structure 
as a two-dimensional resistor network. Precisely, we con- 
sider a square-lattice of N x N resistors, where A'' deter- 
mines the linear size of the network. This lattice lies on 
an insulating substrate at a given temperature Tq, which 
acts as a thermal bath. Each resistor can be in two dif- 
ferent states*^'^^'*^: (i) regular, corresponding to a resis- 
tance r„ ~ ro[l+Q!(r„— To)], and (ii) broken, correspond- 
ing to an effectively "infinite" resistance, rop = 10^r„ 
(resistors in this state will be called defects). In the pre- 
vious expression of r„, a is the temperature coefficient of 
the resistance and T„ the local temperature. This latter 
is determined by Joule heating effects*^'^*'*^ and thermal 
exchanges between neighbor resistors^^: 

T„ = To + A[rnil + {S/ANneig) - 

where, i,,, is the current flowing in the nth resistor and 
Nneig the number of nearest neighbors over which the 
summation is performed. The parameter A represents 



the thermal resistance of each resistor and sets the im- 
portance of Joule heating effects. By taking the above 
expression for T„ we are assuming an instantaneous 
thermalization*^ of each resistor at the value T„. The 
external bias consists of a constant current I applied 
through perfectly conducting bars at the left and right 
sides of the network. In the initial state of the network 
(no external bias) all the resistors are identical: r„ =ro- 
Wc assume that two competing biased processes act 
to determine the evolution of the network^ ^ . These 
two processes consist of stochastic transitions between 
the two possible states of each resistor and they are taken 
to occur through thermal activation, with probabilities^^: 
Won = expl-Eo/kBTn] and Wru = exp[-Ei{/kBTn], 
characterized by the two activation energies, Ejj and Er 
{ks being the Boltzmann constant). The time evolution 
of the network is obtained by Monte Carlo simulations 
which update the network resistance after breaking and 
recovery processes, according to an iterative procedure 
described in details in Rcf The sequence of successive 
configurations provides a resistance signal, R{t), after an 
appropriate calibration of the time scale. Then, depend- 
ing on the stress conditions (/ and Tq) and on the network 
parameters (size, activation energies and other parame- 
ters related to the material like ro and a), the network 
either reaches a steady state or undergoes an irreversible 
electrical failure^*'^^. This latter possibility is associated 
with the achievement of the percolation threshold, Pc, 
for the fraction of broken resistors^^. Therefore, for a 
given network at a given temperature, a threshold cur- 
rent value, Ib, exists above which electrical breakdown 
occurs^*. For values of the current below this thresh- 
old, the steady state of the network is characterized by 
fiuctuations of the fraction of broken resistors, 8p. and of 
the resistance, i5it!, around their respective average values 
<p> and < R>. In particular, we underline that in the 
vanishing current limit (random percolation)^^, the ratio 
{Ed ~ ER)/kBTo determines the average fraction of de- 
fects and thus the level of disorder inside the network. In 
the following we analyze the results of simulations per- 
formed by considering networks of different sizes, with 
different levels of disorder and stressed by different cur- 
rents at room temperature, To = 300 (K). In all the cases 
we take Ed = 0.170 (eV), ro = 1 (0), a = 10'^ (R-^), 
A = 5x 10^ (K/W) (these values are chosen as physically 
reasonable). The values of N range between 100 150, 
while Er between 0.026 ^ 0.164 (eV). 



III. RESULTS 

The resistance of a 100 x 100 network is reported as 
a function of the time in Fig. 1. The different curves, 
starting from the bottom, are obtained for increasing val- 
ues of the external current. In particular, all the evolu- 
tions in Fig. 1 are obtained by taking the activation 
energy of the recovery process equal to Er = 0.103 eV, a 
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value which leads to a network with an intermediate level 
of disorder. The first three curves are associated with 
steady states of the network. More precisely, the lowest 
curve corresponds to the linear response regime and it 
is obtained for / = 0.0013A < Iq, where Iq is the cur- 
rent value associated with the onset of the nonlinearity of 
the I-V characteristic^^. The second curve corresponds 
to the nonlinear regime and it is obtained for I = 0.70 
A. The third corresponds to the threshold for electrical 
breakdown which occurs for I > Ib = 0.95(A). By con- 
trast, the highest curve displays the resistance evolution 
of a network undergoing electrical breakdown and it is 
obtained for a current / = 1.05(A). Overall, Fig. 1 il- 
lustrates qualitatively two important features of the elec- 
trical response of a conducting film. First, it evidences 
that the linear regime occurring for 7 < Jq is followed 
by a nonlinear regime where the average resistance in- 
creases significantly at increasing current^^'^''. Second, 
it shows that the amplitude of the resistance fluctuations 
increases strongly with the external bias. This amplifica- 
tion of the fluctuations becomes particularly important 
when the current reaches the threshold for breakdown 
(third curve from the bottom in Fig. 1). A detailed anal- 
ysis of the behavior of the average resistance and of the 
relative variance of resistance fluctuations as a function 
of the current can be found in Refs.^^'^^ 
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FIG. 1. Resistance evolutions of a 100 x 100 network 
stressed by increasing current values. Starting from the bot- 
tom: I = 0.0013 A (linear regime), I = 0.70 A (nonlinear 
regime), I = 0.95 A (threshold current), I = 1.05 A (break- 
down). The recovery energy is Er = 0.103 eV and corre- 
sponds to an intermediate level of disorder. 

The distributions of the resistance fluctuations, 6R, 
for the three steady state signals in Fig. 1 are reported 
in Fig. 2. Precisely, by denoting with $ the PDF of 
the SR distribution and with a the root mean square 
deviation from the average resistance value, we plot in 
Fig. 2 on a lin-log scale the product C7$ as a function 
of (< i? > —R)/a. This normalized representation, by 
making the distribution independent of its first and sec- 
ond moments, is particularly convenient to explore the 
functional form of any distribution^'^. Here the PDFs 



have been calculated by considering time series contain- 
ing about 1.2 X 10® resistance values. In this figure, the 
data represented by crosses, small circles and up trian- 
gles are obtained respectively for / = 0.0013, 0.70, 0.95 
A (i.e. linear regime, nonlinear regime and threshold 
for breakdown). For comparison, in this figure we also 
report with a dashed curve the Gaussian distribution 
(which in this normalized representation has zero mean 
and unit variance) and, with a continuous curve, the BHP 
distribution^''*. By defining the normalized variable, 
y = {R- < R >)/a, the normalized PDF, U{y) = a^{y) 
and X = b{y — s), the BHP distribution has the following 
functional form^: 



(2) 



where a = 7r/2, b = 0.936 ± 0.002, s = 0.374 ± 0.001 
and if = 2.15 ±0.01"*. This expression can be considered 
as a generalization of the Gumbel distribution, which is 
often associated with the occurrence of rare events. Fig- 
ure 2 shows a remarkable deviation from the Gaussian 
behavior when the network is stressed by the threshold 
current Is- We have found that in this case the PDF 
is well fitted by the Eq. (2), once the parameters a, b, 
s and K are taken as fitting parameters'^ In any case, 
we emphasize the fact that, at least for a network with 
intermediate level of disorder, such as that considered in 
Figs. 1 and 2, the deviations from Gaussianity remain 
weak even for current values in the nonlinear regime and 
become important only close to the electrical breakdown. 
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FIG. 2. Normalized PDF of resistance fluctuations for a 
network of size 100 x 100 biased by / = 0.0013 A (crosses, 
linear regime), I = 0.70 A (small circles, nonlinear regime), 
/ — 0.95 A (up triangles, threshold current). The recovery 
energy is the same of Fig. 1. The thick solid curve and the 
dashed one correspond to the BHP and Gaussian distribu- 
tions, respectively. 

The effect of the size of the network on the distribu- 
tion of resistance fluctuations, is analyzed in Figs. 3 and 
4. Figure 3 compares the PDF obtained for a current 
in the linear regime and applied to a network 100 x 100 
(the data, reported by crosses, are the same of Fig. 2) 
with the PDF obtained for a current in the same regime 
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applied to a network 125 x 125 (diamonds). The two net- 
works have the same parameters and differ only for the 
size, thus they share the same level of internal disorder. 
Figure 3 shows that in the linear regime the distribution 
of SR is Gaussian for all system sizes. In weak 
non-Gaussian tail appears for very small systems that 
vanishes for systems of larger size. 



10" 



10- 
10-^ 
10-' 
lO" 
10"' 



Linear regime 

— Gaussian 

— BHP 




-7 



-5 



-3-11 3 5 
(<R>-R)/a 

FIG. 3. Normalized PDF of resistance fluctuations for net- 
works of size 100 X 100 (crosses) and 125 x 125 (diamonds). 
The networks are biased by currents in the linear regime. The 
recovery energy is the same of Fig. 1. The thick solid curve 
and the dashed one correspond to the BHP and Gaussian 
distributions, respectively. 
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FIG. 4. Normahzcd PDF of resistance fluctuations for net- 
works of size 100 X 100 (up triangles) and 125 x 125 (down tri- 
angles). The networks are biased by / = 0.95 A and / = 1.15 
A which correspond to their respective threshold for break- 
down. The recovery energy Ea is the same of the previous 
figures and the thick solid and the dashed curves have the 
same meaning of Fig. 3. 

Figure 4 displays a comparison similar to that in Fig. 
3 for two networks of size 100 x 100 (up triangles, same 
data of Fig. 2) and 125x125 (down triangles), stressed by 
currents corresponding to the threshold for breakdown. 
The current values are / = 0.95, 1.15 A, respectively for 



the networks with TV = 100 and = 125. We can see 
that in both cases the PDFs exhibit non-Gaussian tails. 
However, the same figure shows that the non-Gaussianity 
is weaker for the system of larger size. By performing 
a detailed investigation on systems of different sizes we 
have found that this trend is systematic^*^ . Therefore, we 
conclude that, at least when networks with intermediate 
level of disorder are considered, as in the case of Figs. 
1-4, the deviations from Gaussianity are related to the 
finite size of the system. 

Now, before discussing the role of disorder in the break- 
down process, it is convenient to distinguish between the 
intrinsic disorder and the disorder driven by the external 
bias. For given values of Tq and Ejj, the average fraction 
of broken resistors in the vanishing current limit, < p >o, 
is only determined by the recovery energy En^^ . This av- 
erage fraction of defects represents an intrinsic property 
of the system, which characterizes the conducting mate- 
rial and sets the level of disorder inside the network in ab- 
sence of an external bias. We note nate the steady state 
of a network in this vanishing current limit is determined 
by the competition between two random percolations^^. 
On the other hand, for I > Iq, there is also an additional 
component of disorder due to the defects generated by 
the external current^^. Therefore, the average fraction 
of defects < p > becomes dependent on the current and 
its relative variation, [< p > — < p >o]/ < p >o, scales 
as^^ {I/Io)^. Of course, this dependence of < p > on 
/ is important because it actually drives the system to- 
wards the breakdown. However, it is also crucial to point 
out the role of the intrinsic disorder on the breakdown 
process and its effect on the distribution of resistance 
fluctuations. Therefore, in the following we will consider 
networks with different values of the recovery activation 
energy Er. 
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FIG. 5. Normalized PDF of resistance fiuctuations of 
100 X 100 networks at increasing Er values (level of dis- 
order). In all the cases the current value corresponds to 
the respective threshold for breakdown. Precisely, squares: 
Er = 0.026, Ib = 3.5; up triangles: Er = 0.103, Ib = 0.95; 
stars: Er = 0.155, Ib = 0.11, left triangles: Er = 0.164, 
Ib = 0.009 (energies in eV, currents in A). 
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In a previous work^® we have shown that for an ar- 
bitrary value of En satisfing the stabiUty condition^®'^^, 
i.e. providing a steady state of the network, the elec- 
trical breakdown is associated with a first order phase 
transition. This result agrees with the behavior observed 
in electrical breakdown experiments, performed in the 
Joule regime of composites^^'^*. Nevertheless, it can 
be shown^'^ that when En reaches its maximum value 
compatible with a steady state, E^^max, the conductor- 
insulator transition becomes of the second order. This 
change in the order of the transition, when going from 
low to high disordered systems, has been also predicted 
by Andersen et al.^^ 

Therefore, we report in Fig. 5 the PDFs of resistance 
fluctuations of 100 x 100 networks characterized by dif- 
ferent values of Er. The PDFs have been calculated for 
current values corresponding to the breakdown threshold. 
We can see that the non-Gaussianity of the distribution 
at I = Ib increases systematically at increasing values 
of En (i.e. at increasing level of the intrinsic disorder). 
In particular. Fig. 5 shows that when the value of En is 
very close to En.MAx, and thus the system approaches 
the critical point, the PDF achieves (left triangles) the 
BHP form. In these conditions, the PDF should become 
independent of the system size^'*. 




(<R>-R)/o 

FIG. 6. Normalized PDF of resistance fluctuations for net- 
works of increasing size and close to the critical conditions: 
the recovery energy is Er = 0.164 eV and the bias current 
corresponds to the threshold for breakdown. Precisely, left 
triangles refer to a network 100 x 100 biased by 7 = 0.0090 
A, circles to a network 125 x 125 biased by 7 = 0.0011 A and 
right triangles to a network 150 x 150 and 7 = 0.0013 A. 

Indeed, Fig. 6 confirms this behavior. In this figure we 
report the PDFs calculated for the three networks of size 
100 X 100 (left triangles, same data of Fig. 5), 125 x 125 
(circles) and 150 x 150 (right triangles). In all the cases 
the recovery energy is En = O.lGAeV « En,MAX and the 
external current corresponds to the threshold value for 
breakdown {Ib = 0.0090,0.011,0.013 A, respectively). 
The data in Fig. 6 shows that, when En is very close 
to En,MAx and thus the system approaches the critical 



point, the PDFs of resistance fluctuations become inde- 
pendent of the system size and they are well described 
by the BHP distribution. 
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FIG. 7. Power spectral density of resistance fluctuations of 
125 X 125 networks for increasing Er values. Precisely, start- 
ing from the bottom: Er = 0.026,0.103,0.155,0.164 eV. In 
all cases the current corresponds to the respective threshold 
for breakdown (respectively: 7 = 4.4, 1.15, 0.14, 0.011 A). For 
visual reasons the lowest spectrum has been shifted down by 
multipling it by a factor 5 x 10~^. The gray lines represent 
the best-fit with a Lorentzian spectral density. 

Other interesting information about the system can 
be extracted from the spectral analysis of resistance 
fluctuations^'^'^. Hence we report in Fig. 7 the power 
spectral density of resistance fluctuations of 125 x 125 
networks at increasing values of En- AH the spectra are 
calculated for networks stressed by currents correspond- 
ing to the threshold for electrical breakdown. Starting 
from the bottom, the curves arc obtained respectively for 
En = 0.026,0.103,0.155 and 0.164 eV (the same values 
of En considered in Fig. 5). Basically all the spectra dis- 
play a Lorentzian behavior which witnesses an exponen- 
tial decay of the auto-correlation function of resistance 
fluctuations^ (the gray curves in this figure represent the 
best-fit with a Lorentzian distribution). Moreover, we 
can sec that for increasing values of En, i.e. for increas- 
ing level of intrinsic disorder, the corner frequency of 
the spectra is systematically shifted towards lower val- 
ues. This occurs because, at increasing values of < p >, 
the network approaches the percolation threshold, thus 
the correlation length increases as ^ ^ | < p > —pc\^'^ 
(where v is the correlation length exponent)^^ and the 
same occurs for the correlation time of resistance fluctu- 
ations. However, we note that for the two upper spectra 
in Fig. 7 the flt with a Lorentzian curve appears progres- 
sively less satisfactory in the low frequency region. On 
the other hand, these spectra are obtained by taking the 
value of En near or very near to En,MAX and thus corre- 
spond to networks progressively closer to the critical con- 
ditions. This anomalous behavior in the low frequency 
region of the power specral density of resistance fluctu- 
ations in the case of networks close to criticality, could 
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even suggest the emergence of a 1// behavior. Longer 
time series are required to point out the behavior of the 
spectra in this frequency region. 



IV. CONCLUSIONS 

We have studied the distribution of the resistance fluc- 
tuations of conducting thin films with different levels of 
internal disorder. The study has been performed by de- 
scribing the film as a resistor network in a steady state 
determined by the competition of two biased stochastic 
processes, according to the SBRN modeP°~^^. We have 
considered systems of different sizes and under different 
stress conditions, from the linear response regime up to 
the threshold for electrical breakdown. A remarkable 
non-Gaussianity of the fluctuation distribution is found 
near breakdown. This non-Gaussianity becomes more ev- 
ident at increasing the degree of disorder of the network. 
As a general trend, these deviations from Gaussianity 
are related to the finite size of the system and they are 
found to vanish in the large size limit. However, near the 
critical point of the conductor-insulator transition, the 
non-Gaussianity is found to persist in the large size limit 
and is well fitted by the universal Bramwell-Holdsworth- 
Pinton distribution^'^. Furthermore, in the last case, the 
non-Gaussianity seems to be associated with an anoma- 
lous behavior in the low frequency region of the specral 
density of resistance fluctuations. 
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